Figures for scRNAseq data
Anne Hoffrichter
2021/01/22
1 Goal
Preparation of figures for publication
2 Approach
library(Seurat)
library(viridis)
library(org.Hs.eg.db)
library(AnnotationDbi)
library(tidyverse)
library(purrr)
library(patchwork)
library(grid)
library(clusterProfiler)
set.seed(42)boost <- readRDS("../data/boost_data_seurat.rds")
cds <- readRDS("../data/boost_data_cds.rds")3 Figures
cols.cluster <- hcl.colors(6, "Spectral") %>% rev()
cols.cluster <- cols.cluster[c(1, 3, 4, 5, 6, 2)]3.1 Figure 1 F: Dimensional reduction plot in UMAP space, cells colored by cluster membership, or pseudotime
DimPlot(boost, cols = cols.cluster, pt.size = 2) + theme(axis.text = element_text(size = 25),
axis.title = element_text(size = 30), plot.title = element_blank(), legend.text = element_text(size = 20)) +
guides(colour = guide_legend(override.aes = list(size = 8))) + scale_y_continuous(breaks = c(-5,
0, 5))
FeaturePlot(boost, features = "pseudotime", cols = plasma(length(boost$pseudotime)),
pt.size = 2) + theme(axis.text = element_text(size = 25), axis.title = element_text(size = 30),
plot.title = element_blank(), legend.text = element_text(size = 20), legend.key.size = unit(0.8,
"cm")) + scale_y_continuous(breaks = c(-5, 0, 5))## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.2 Figure 1 G: Heatmap of specific marker genes. Cells are ordered by pseudotime, with an additional colorbar for the respective cluster.
genes <- c("MKI67", "SOX2", "HES1", "HES5", "NES", "PAX6", "TUBB3", "STMN1", "FAT3",
"DCX", "RBFOX3", "MAPT", "SNAP25", "CAMK2A", "CAMK4", "NRXN1", "NCAM1", "DLG4",
"SYP", "SYN1", "SYT1", "ANK3", "SPTBN4", "SCN2A", "GRIN1", "GRIN2B", "GRIA1",
"GRIA2", "GABRA3", "GABRB2", "ARC", "NPAS4", "FOXG1")DoMultiBarHeatmap2 <- function(object, features = NULL, cells = NULL, group.by = "ident",
additional.group.by = NULL, additional.group.sort.by = NULL, cols.use = NULL,
group.bar = TRUE, disp.min = -2.5, disp.max = NULL, slot = "scale.data", assay = NULL,
label = TRUE, size = 5.5, hjust = 0, angle = 45, raster = TRUE, draw.lines = TRUE,
lines.width = NULL, group.bar.height = 0.02, combine = TRUE, group.name = NULL,
additional.group.name = NULL) {
cells <- cells %||% colnames(x = object)
if (is.numeric(x = cells)) {
cells <- colnames(x = object)[cells]
}
assay <- assay %||% DefaultAssay(object = object)
DefaultAssay(object = object) <- assay
features <- features %||% VariableFeatures(object = object)
features <- rev(x = unique(x = features))
disp.max <- disp.max %||% ifelse(test = slot == "scale.data", yes = 2.5, no = 6)
possible.features <- rownames(x = GetAssayData(object = object, slot = slot))
if (any(!features %in% possible.features)) {
bad.features <- features[!features %in% possible.features]
features <- features[features %in% possible.features]
if (length(x = features) == 0) {
stop("No requested features found in the ", slot, " slot for the ", assay,
" assay.")
}
warning("The following features were omitted as they were not found in the ",
slot, " slot for the ", assay, " assay: ", paste(bad.features, collapse = ", "))
}
if (!is.null(additional.group.sort.by)) {
if (any(!additional.group.sort.by %in% additional.group.by)) {
bad.sorts <- additional.group.sort.by[!additional.group.sort.by %in%
additional.group.by]
additional.group.sort.by <- additional.group.sort.by[additional.group.sort.by %in%
additional.group.by]
if (length(x = bad.sorts) > 0) {
warning("The following additional sorts were omitted as they were not a subset of additional.group.by : ",
paste(bad.sorts, collapse = ", "))
}
}
}
data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, slot = slot)[features,
cells, drop = FALSE])))
object <- suppressMessages(expr = StashIdent(object = object, save.name = "ident"))
group.by <- group.by %||% "ident"
groups.use <- object[[c(group.by, additional.group.by[!additional.group.by %in%
group.by])]][cells, , drop = FALSE]
plots <- list()
for (i in group.by) {
data.group <- data
if (!is_null(additional.group.by)) {
additional.group.use <- additional.group.by[additional.group.by != i]
if (!is_null(additional.group.sort.by)) {
additional.sort.use = additional.group.sort.by[additional.group.sort.by !=
i]
} else {
additional.sort.use = NULL
}
} else {
additional.group.use = NULL
additional.sort.use = NULL
}
group.use <- groups.use[, c(i, additional.group.use), drop = FALSE]
for (colname in colnames(group.use)) {
if (!is.factor(x = group.use[[colname]])) {
group.use[[colname]] <- factor(x = group.use[[colname]])
}
}
if (draw.lines) {
lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 0.0025)
placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) *
lines.width), FUN = function(x) {
return(Seurat:::RandomName(length = 20))
})
placeholder.groups <- data.frame(rep(x = levels(x = group.use[[i]]),
times = lines.width))
group.levels <- list()
group.levels[[i]] = levels(x = group.use[[i]])
for (j in additional.group.use) {
group.levels[[j]] <- levels(x = group.use[[j]])
placeholder.groups[[j]] = NA
}
colnames(placeholder.groups) <- colnames(group.use)
rownames(placeholder.groups) <- placeholder.cells
group.use <- sapply(group.use, as.vector)
rownames(x = group.use) <- cells
group.use <- rbind(group.use, placeholder.groups)
for (j in names(group.levels)) {
group.use[[j]] <- factor(x = group.use[[j]], levels = group.levels[[j]])
}
na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells),
ncol = ncol(x = data.group), dimnames = list(placeholder.cells, colnames(x = data.group)))
data.group <- rbind(data.group, na.data.group)
}
order_expr <- paste0("order(", paste(c(i, additional.sort.use), collapse = ","),
")")
group.use = with(group.use, group.use[eval(parse(text = order_expr)), , drop = F])
plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, disp.min = disp.min,
disp.max = disp.max, feature.order = features, cell.order = rownames(x = group.use),
group.by = group.use[[i]])
if (group.bar) {
pbuild <- ggplot_build(plot = plot)
group.use2 <- group.use
cols <- list()
na.group <- Seurat:::RandomName(length = 20)
for (colname in rev(x = colnames(group.use2))) {
if (colname == i) {
colid = group.name
} else {
colid = additional.group.name
}
# Default
cols[[colname]] <- c((scales::hue_pal())(length(x = levels(x = group.use[[colname]]))))
# Overwrite if better value is provided
if (!is_null(cols.use[[colname]])) {
req_length = length(x = levels(group.use))
if (length(cols.use[[colname]]) < req_length) {
warning("Cannot use provided colors for ", colname, " since there aren't enough colors.")
} else {
if (!is_null(names(cols.use[[colname]]))) {
if (all(levels(group.use[[colname]]) %in% names(cols.use[[colname]]))) {
cols[[colname]] <- as.vector(cols.use[[colname]][levels(group.use[[colname]])])
} else {
warning("Cannot use provided colors for ", colname, " since all levels (",
paste(levels(group.use[[colname]]), collapse = ","), ") are not represented.")
}
} else {
cols[[colname]] <- as.vector(cols.use[[colname]])[c(1:length(x = levels(x = group.use[[colname]])))]
}
}
}
# Add white if there's lines
if (draw.lines) {
levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]),
na.group)
group.use2[placeholder.cells, colname] <- na.group
cols[[colname]] <- c(cols[[colname]], "#FFFFFF")
}
names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])
y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range *
0.015
y.max <- y.pos + group.bar.height * y.range
pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1],
y.max)
plot <- suppressMessages(plot + annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]),
xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + annotation_custom(grob = grid::textGrob(label = colid,
hjust = 0, gp = gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)),
ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) + coord_cartesian(ylim = c(0,
y.max), clip = "off"))
if ((colname == i) && label) {
x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
x.divs <- pbuild$layout$panel_params[[1]]$x.major %||% pbuild$layout$panel_params[[1]]$x$break_positions()
group.use$x <- x.divs
label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]],
FUN = median) * x.max
label.x.pos <- data.frame(group = names(x = label.x.pos), label.x.pos)
plot <- plot + geom_text(stat = "identity", data = label.x.pos,
aes_string(label = "group", x = "label.x.pos"), y = y.max + y.max *
0.03 * 0.5, angle = angle, hjust = hjust, size = size)
plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, y.max +
y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) *
size), clip = "off"))
}
}
}
plot <- plot + theme(line = element_blank())
plots[[i]] <- plot
}
if (combine) {
plots <- CombinePlots(plots = plots)
}
return(plots)
}boost_subset <- subset(boost, new_order %in% c("1", "2", "3", "4", "5"))
cols.use <- list(pseudotime = plasma(length(boost_subset$pseudotime)), ident = cols.cluster)DoMultiBarHeatmap2(boost_subset, features = genes, group.by = "pseudotime", additional.group.by = "ident",
draw.lines = FALSE, label = FALSE, cols.use = cols.use, raster = FALSE) + scale_fill_viridis_c(option = "viridis") +
NoLegend() + theme(axis.text = element_text(size = 14, colour = "black")) + scale_y_discrete(position = "left") +
geom_hline(yintercept = c(3.5, 9.5, 12.5, 16.5, 23.5, 27.5), color = "white",
size = 1)## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.# DoHeatmap(boost_subset, features = genes, group.by='ident', draw.lines = FALSE,
# label=FALSE, raster = FALSE)+ scale_fill_viridis_c(option='viridis')+
# theme(legend.text = element_text(size = 20), legend.key.size = unit(0.8, 'cm'),
# legend.text.align = 1, legend.title = element_text(size=30))3.3 Figure 1 H: NMI scores per cluster of neuronal cells.
sub <- subset(boost, idents = c("2", "3", "4", "5"))a <- VlnPlot(sub, features = "NMIoverall", pt.size = 0.1, cols = cols.cluster[2:5]) +
NoLegend() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ggtitle("NMI") + ylab("")
b <- VlnPlot(sub, features = "NMIdiscriminable", pt.size = 0.1, cols = cols.cluster[2:5]) +
NoLegend() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ggtitle("Discriminating NMI") +
ylab("")
c <- VlnPlot(sub, features = "NMIactivity", pt.size = 0.1, cols = cols.cluster[2:5]) +
NoLegend() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ggtitle("Functionality index") +
ylab("")# pw <- (a/b/c)&theme(axis.text.x = element_text(angle = 0, hjust =
# 0.5),axis.text=element_text(size=20),
# plot.title=element_text(size=25))&xlab('')# pw&theme(axis.text.x = element_text(angle = 0, hjust =
# 0.5),axis.text=element_text(size=20),
# plot.title=element_text(size=22))&xlab('')&ylab('')b + scale_y_continuous(breaks = c(0.49, 0.52, 0.55))## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
c
3.4 Figure 2 D: Violin Plots of CASP9, CASP3, CASP7
genes <- c("CASP7", "CASP9", "CASP3")p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster, combine = F) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)3.4.1 CASP9
p[[2]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0,
hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_blank()) +
xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.4.2 CASP3
p[[3]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0",
"3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.4.3 CASP7
p[[1]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()
3.5 Figure 2 I: Violin Plots BAX, BCL2, BAK1
genes <- c("BAX", "BCL2", "BAK1")p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))3.5.1 BAX
p[[1]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.5.2 BCL2
p[[2]] + scale_y_continuous(breaks = c(0, 0.4, 0.8, 1.2)) + theme(axis.text.x = element_text(angle = 0,
hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) +
xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.5.3 BAK1
p[[3]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0), labels=c('0.0', '1.0', '2.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.6 Figure 2 M: Violin Plot CYCS
genes <- "CYCS"p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))p[[1]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0",
"3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.7 Figure 3 C: Violin Plots BIRC2, BIRC5, AREL1, BIRC4 (XIAP)
genes <- c("BIRC2", "XIAP", "BIRC5", "AREL1")p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster, combine = F) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))3.7.1 BIRC2
p[[1]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.7.2 BIRC5
p[[3]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.7.3 AREL1
p[[4]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.7.4 XIAP/BIRC4
p[[2]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0), labels=c('0.0', '1.0', '2.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.8 Figure 3 F: Violin Plots HTRA2, SIVA1, SIAH1
genes <- c("HTRA2", "SIVA1", "SIAH1", "USP11")p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))3.8.1 HTRA2
p[[1]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.8.2 SIVA1
p[[2]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.8.3 SIAH1
p[[3]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.9 Figure 3 G: Violin Plot USP11
p[[4]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0",
"3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.10 Supplementary Figure 1 H: Pseudotime trajectories calculated on the per cell expression values for specific genes.
my_plot_genes <- function(cds_subset, min_expr = NULL, cell_size = 0.75, nrow = NULL,
ncol = 1, panel_order = NULL, color_cells_by = "pseudotime", trend_formula = "~ splines::ns(pseudotime, df=3)",
label_by_short_name = TRUE, vertical_jitter = NULL, horizontal_jitter = NULL,
line.size = NULL) {
assertthat::assert_that(methods::is(cds_subset, "cell_data_set"))
tryCatch({
pseudotime(cds_subset)
}, error = function(x) {
stop(paste("No pseudotime calculated. Must call order_cells first."))
})
colData(cds_subset)$pseudotime <- pseudotime(cds_subset)
if (!is.null(min_expr)) {
assertthat::assert_that(assertthat::is.number(min_expr))
}
assertthat::assert_that(assertthat::is.number(cell_size))
if (!is.null(nrow)) {
assertthat::assert_that(assertthat::is.count(nrow))
}
assertthat::assert_that(assertthat::is.count(ncol))
assertthat::assert_that(is.logical(label_by_short_name))
if (label_by_short_name) {
assertthat::assert_that("gene_short_name" %in% names(rowData(cds_subset)),
msg = paste("When label_by_short_name = TRUE,", "rowData must have a column of gene",
"names called gene_short_name."))
}
assertthat::assert_that(color_cells_by %in% c("cluster", "partition") | color_cells_by %in%
names(colData(cds_subset)), msg = paste("color_cells_by must be a column in the",
"colData table."))
if (!is.null(panel_order)) {
if (label_by_short_name) {
assertthat::assert_that(all(panel_order %in% rowData(cds_subset)$gene_short_name))
} else {
assertthat::assert_that(all(panel_order %in% row.names(rowData(cds_subset))))
}
}
assertthat::assert_that(nrow(rowData(cds_subset)) <= 100, msg = paste("cds_subset has more than 100 genes -",
"pass only the subset of the CDS to be", "plotted."))
assertthat::assert_that(methods::is(cds_subset, "cell_data_set"))
assertthat::assert_that("pseudotime" %in% names(colData(cds_subset)), msg = paste("pseudotime must be a column in",
"colData. Please run order_cells", "before running", "plot_genes_in_pseudotime."))
if (!is.null(min_expr)) {
assertthat::assert_that(assertthat::is.number(min_expr))
}
assertthat::assert_that(assertthat::is.number(cell_size))
assertthat::assert_that(!is.null(size_factors(cds_subset)))
if (!is.null(nrow)) {
assertthat::assert_that(assertthat::is.count(nrow))
}
assertthat::assert_that(assertthat::is.count(ncol))
assertthat::assert_that(is.logical(label_by_short_name))
if (label_by_short_name) {
assertthat::assert_that("gene_short_name" %in% names(rowData(cds_subset)),
msg = paste("When label_by_short_name = TRUE,", "rowData must have a column of gene",
"names called gene_short_name."))
}
assertthat::assert_that(color_cells_by %in% c("cluster", "partition") | color_cells_by %in%
names(colData(cds_subset)), msg = paste("color_cells_by must be a column in the",
"colData table."))
if (!is.null(panel_order)) {
if (label_by_short_name) {
assertthat::assert_that(all(panel_order %in% rowData(cds_subset)$gene_short_name))
} else {
assertthat::assert_that(all(panel_order %in% row.names(rowData(cds_subset))))
}
}
assertthat::assert_that(nrow(rowData(cds_subset)) <= 100, msg = paste("cds_subset has more than 100 genes -",
"pass only the subset of the CDS to be", "plotted."))
f_id <- NA
Cell <- NA
cds_subset = cds_subset[, is.finite(colData(cds_subset)$pseudotime)]
cds_exprs <- SingleCellExperiment::counts(cds_subset)
cds_exprs <- Matrix::t(Matrix::t(cds_exprs)/size_factors(cds_subset))
cds_exprs <- reshape2::melt(round(as.matrix(cds_exprs)))
if (is.null(min_expr)) {
min_expr <- 0
}
colnames(cds_exprs) <- c("f_id", "Cell", "expression")
cds_colData <- colData(cds_subset)
cds_rowData <- rowData(cds_subset)
cds_exprs <- as.data.frame(merge(cds_exprs, cds_rowData, by.x = "f_id", by.y = "row.names"))
cds_exprs <- as.data.frame(merge(cds_exprs, cds_colData, by.x = "Cell", by.y = "row.names"))
cds_exprs$adjusted_expression <- cds_exprs$expression
if (label_by_short_name == TRUE) {
if (is.null(cds_exprs$gene_short_name) == FALSE) {
cds_exprs$feature_label <- as.character(cds_exprs$gene_short_name)
cds_exprs$feature_label[is.na(cds_exprs$feature_label)] <- cds_exprs$f_id
} else {
cds_exprs$feature_label <- cds_exprs$f_id
}
} else {
cds_exprs$feature_label <- cds_exprs$f_id
}
cds_exprs$f_id <- as.character(cds_exprs$f_id)
cds_exprs$feature_label <- factor(cds_exprs$feature_label)
new_data <- data.frame(pseudotime = colData(cds_subset)$pseudotime)
model_tbl = fit_models(cds_subset, model_formula_str = trend_formula)
model_expectation <- model_predictions(model_tbl, new_data = colData(cds_subset))
colnames(model_expectation) <- colnames(cds_subset)
expectation <- plyr::ddply(cds_exprs, plyr::.(f_id, Cell), function(x) {
data.frame(expectation = model_expectation[x$f_id, x$Cell])
})
cds_exprs <- merge(cds_exprs, expectation)
cds_exprs$expression[cds_exprs$expression < min_expr] <- min_expr
cds_exprs$expectation[cds_exprs$expectation < min_expr] <- min_expr
if (!is.null(panel_order)) {
cds_exprs$feature_label <- factor(cds_exprs$feature_label, levels = panel_order)
}
q <- ggplot(aes(pseudotime, expression), data = cds_exprs)
if (!is.null(color_cells_by)) {
q <- q + geom_point(aes_string(color = color_cells_by), size = I(cell_size),
position = position_jitter(horizontal_jitter, vertical_jitter))
if (class(colData(cds_subset)[, color_cells_by]) == "numeric") {
q <- q + viridis::scale_color_viridis(option = "C")
}
} else {
q <- q + geom_point(size = I(cell_size), position = position_jitter(horizontal_jitter,
vertical_jitter))
}
q <- q + geom_line(aes(x = pseudotime, y = expectation), data = cds_exprs, size = line.size)
q <- q + scale_y_log10() + facet_wrap(~feature_label, nrow = nrow, ncol = ncol,
scales = "free_y")
if (min_expr < 1) {
q <- q + expand_limits(y = c(min_expr, 1))
}
q <- q + ylab("Expression")
q <- q + xlab("pseudotime")
# q <- q + monocle_theme_opts()
q
}genes <- c("SOX2", "NES", "TUBB3", "DCX", "SNAP25", "SYN1")
cds_subset <- cds[, SummarizedExperiment::colData(cds)$ident %in% c("1", "2", "3",
"4", "5")]## Loading required package: monocle3## Loading required package: SingleCellExperiment## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following object is masked from 'package:dplyr':
##
## count## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## The following object is masked from 'package:Biobase':
##
## rowMedians## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb##
## Attaching package: 'SummarizedExperiment'## The following object is masked from 'package:Seurat':
##
## Assays##
## Attaching package: 'monocle3'## The following objects are masked from 'package:Biobase':
##
## exprs, fData, fData<-, pData, pData<-lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "SOX2", ]
p <- list()
p[[1]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4,
line.size = 2) + theme_grey(base_size = 40, ) + NoLegend() + scale_color_manual(values = c("#584B9F",
"#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(),
axis.text.x = element_blank(), strip.background = element_rect(fill = "white",
colour = NA)) + scale_y_log10(breaks = c(0.01, 0.1, 1, 10), labels = c("0.01",
"0.1", "1", "10"))## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "NES", ]
p[[2]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4,
line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F",
"#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title.x = element_blank(),
axis.text.x = element_blank(), strip.background = element_rect(fill = "white",
colour = NA))lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "TUBB3", ]
p[[3]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4,
line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F",
"#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(),
strip.background = element_rect(fill = "white", colour = NA))lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "DCX", ]
p[[4]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4,
line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F",
"#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(),
axis.text.x = element_blank(), strip.background = element_rect(fill = "white",
colour = NA))lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "SNAP25", ]
p[[5]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4,
line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F",
"#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(),
axis.text.x = element_blank(), strip.background = element_rect(fill = "white",
colour = NA))lineage_cds <- cds_subset[rowData(cds_subset)$gene_short_name %in% "SYN1", ]
p[[6]] <- my_plot_genes(lineage_cds, color_cells_by = "new_order", cell_size = 4,
line.size = 2) + theme_grey(base_size = 40) + NoLegend() + scale_color_manual(values = c("#584B9F",
"#BAEEAE", "#FCDE85", "#ED820A", "#A71B4B")) + theme(axis.title = element_blank(),
strip.background = element_rect(fill = "white", colour = NA))wrap_plots(p[1:3], ncol = 1)
wrap_plots(p[4:6], ncol = 1)
3.11 Supplementary Figure 2 B: Violin Plots CASP8, CASP10, CASP2, CASP6
genes <- c("CASP8", "CASP10", "CASP2", "CASP6")p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster, combine = F) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)3.11.1 CASP8
p[[1]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()
3.11.2 CASP10
p[[2]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()
3.11.3 CASP2
p[[3]] + # scale_y_continuous(breaks=c(0.0, 0.3, 0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()
3.11.4 CASP6
p[[4]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()
3.12 Supplementary Figure 2 C: Violin Plot ACIN1
genes <- c("ACIN1")p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster, combine = FALSE)p[[1]] + scale_y_continuous(breaks = c(0, 1, 2), labels = c("0.0", "1.0", "2.0")) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_blank()) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.13 Supplementary Figure 3 E: Violin Plot XAF1
genes <- c("XAF1")p <- VlnPlot(boost_subset, features = genes, group.by = "new_order", pt.size = 0.1,
cols = cols.cluster) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))p[[1]] + scale_y_continuous(breaks = c(0, 0.4, 0.8)) + theme(axis.text.x = element_text(angle = 0,
hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_blank()) +
xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.14 Supplementary Figure 4 C: Glut/Gaba split violin Plots CASP3, CASP9, CASP6, BIRC2, XIAP, AREL1, BCL2, BAK1, BAX, USP11, SIVA1, CYCS
rm(list = setdiff(ls(), c("boost", "cols.cluster")))SLC17A7 <- WhichCells(boost, expression = SLC17A7 > 0)
SLC17A6 <- WhichCells(boost, expression = SLC17A6 > 0)
TBR1 <- WhichCells(boost, expression = TBR1 > 0)
Glut <- union(SLC17A6, SLC17A7)
Glut <- union(Glut, TBR1)GAD1 <- WhichCells(boost, expression = GAD1 > 0)
GAD2 <- WhichCells(boost, expression = GAD2 > 0)
GAB <- union(GAD1, GAD2)both_exp <- FetchData(boost, vars = c("SLC17A7", "SLC17A6", "TBR1", "GAD1", "GAD2"),
cells = intersect(Glut, GAB))
both_exp <- data.frame(Glut = rowMeans(both_exp[, 1:3]), Gaba = rowMeans(both_exp[,
4:5]))
both_exp$cell_id <- rownames(both_exp)
both_exp <- both_exp %>% mutate(type = if_else(Glut > Gaba, "Glut", "Gaba"))
both_GAB <- filter(both_exp, type == "Gaba")
both_Glut <- filter(both_exp, type == "Glut")GAB <- setdiff(GAB, Glut)
GAB <- c(GAB, both_GAB$cell_id)
Glut <- setdiff(Glut, GAB)
Glut <- c(Glut, both_Glut$cell_id)boost[["old.ident"]] <- Idents(object = boost)
boost <- StashIdent(object = boost, save.name = "old.ident")## With Seurat 3.X, stashing identity classes can be accomplished with the following:
## boost[["old.ident"]] <- Idents(object = boost)Idents(object = boost, cells = GAB) <- "Gaba"
Idents(object = boost, cells = Glut) <- "Glut"
Rest <- WhichCells(boost, invert = TRUE, cells = c(GAB, Glut))
Idents(object = boost, cells = Rest) <- "rest"boost_subset <- subset(boost, new_order %in% c("3", "4", "5"))subset <- SubsetData(boost_subset, ident.use = c("Gaba", "Glut"))genes <- c("CASP9", "CASP3", "CASP6")p <- VlnPlot(subset, features = genes, pt.size = 0, cols = cols.cluster, combine = F) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)3.14.1 CASP3
p[[2]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.14.2 CASP9
p[[1]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0,
hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) +
xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.14.3 CASP6
p[[3]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
genes <- c("BIRC2", "XIAP", "AREL1", "SIVA1", "USP11", "BAX", "BCL2", "BAK1", "CYCS")p <- VlnPlot(subset, features = genes, pt.size = 0, cols = cols.cluster, combine = F) #+NoLegend()+xlab('')#+theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
# patchwork <- wrap_plots(p, byrow = T)3.14.4 BIRC2
p[[1]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0,
hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) +
xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.14.5 XIAP
p[[2]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.14.6 AREL1
p[[3]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0",
"3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.14.7 BCL2
p[[7]] + # scale_y_continuous(breaks=c(0.0, 0.3,0.6, 0.9))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.14.8 BAK1
p[[8]] + scale_y_continuous(breaks = c(0, 1, 2, 3), labels = c("0.0", "1.0", "2.0",
"3.0")) + theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.14.9 BAX
p[[6]] + scale_y_continuous(breaks = c(0, 0.3, 0.6, 0.9)) + theme(axis.text.x = element_text(angle = 0,
hjust = 0.5), axis.text = element_text(size = 20), plot.title = element_text(size = 18)) +
xlab("") + ylab("") + NoLegend()## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
3.14.10 USP11
p[[5]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.14.11 SIVA1
p[[4]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
3.14.12 CYCS
p[[9]] + # scale_y_continuous(breaks=c(0.0, 1.0, 2.0, 3.0), labels=c('0.0', '1.0', '2.0',
# '3.0'))+
theme(axis.text.x = element_text(angle = 0, hjust = 0.5), axis.text = element_text(size = 20),
plot.title = element_text(size = 18)) + xlab("") + ylab("") + NoLegend()
4 SessionInfo
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] monocle3_0.2.3.0 SingleCellExperiment_1.12.0
## [3] SummarizedExperiment_1.20.0 GenomicRanges_1.42.0
## [5] GenomeInfoDb_1.26.1 MatrixGenerics_1.2.0
## [7] matrixStats_0.57.0 clusterProfiler_3.18.0
## [9] patchwork_1.1.0 forcats_0.5.0
## [11] stringr_1.4.0 dplyr_1.0.2
## [13] purrr_0.3.4 readr_1.4.0
## [15] tidyr_1.1.2 tibble_3.0.4
## [17] ggplot2_3.3.2 tidyverse_1.3.0
## [19] org.Hs.eg.db_3.12.0 AnnotationDbi_1.52.0
## [21] IRanges_2.24.0 S4Vectors_0.28.0
## [23] Biobase_2.50.0 BiocGenerics_0.36.0
## [25] viridis_0.5.1 viridisLite_0.3.0
## [27] Seurat_3.2.2
##
## loaded via a namespace (and not attached):
## [1] reticulate_1.18 tidyselect_1.1.0 RSQLite_2.2.1
## [4] htmlwidgets_1.5.2 BiocParallel_1.24.1 Rtsne_0.15
## [7] scatterpie_0.1.5 munsell_0.5.0 codetools_0.2-18
## [10] ica_1.0-2 future_1.20.1 miniUI_0.1.1.1
## [13] withr_2.3.0 colorspace_2.0-0 GOSemSim_2.16.1
## [16] knitr_1.30 rstudioapi_0.13 ROCR_1.0-11
## [19] tensor_1.5 DOSE_3.16.0 listenv_0.8.0
## [22] labeling_0.4.2 GenomeInfoDbData_1.2.4 polyclip_1.10-0
## [25] bit64_4.0.5 farver_2.0.3 downloader_0.4
## [28] parallelly_1.21.0 vctrs_0.3.5 generics_0.1.0
## [31] xfun_0.31 R6_2.5.0 graphlayouts_0.7.1
## [34] rsvd_1.0.3 DelayedArray_0.16.0 bitops_1.0-6
## [37] spatstat.utils_1.17-0 fgsea_1.16.0 assertthat_0.2.1
## [40] promises_1.1.1 scales_1.1.1 ggraph_2.0.4
## [43] enrichplot_1.10.1 gtable_0.3.0 globals_0.14.0
## [46] goftest_1.2-2 tidygraph_1.2.0 rlang_1.0.3
## [49] splines_4.0.3 lazyeval_0.2.2 broom_0.7.2
## [52] BiocManager_1.30.10 yaml_2.2.1 reshape2_1.4.4
## [55] abind_1.4-5 modelr_0.1.8 backports_1.2.0
## [58] httpuv_1.5.4 qvalue_2.22.0 tools_4.0.3
## [61] bookdown_0.21 ellipsis_0.3.1 jquerylib_0.1.4
## [64] RColorBrewer_1.1-2 ggridges_0.5.2 Rcpp_1.0.5
## [67] plyr_1.8.6 zlibbioc_1.36.0 RCurl_1.98-1.2
## [70] rpart_4.1-15 deldir_0.2-3 pbapply_1.4-3
## [73] cowplot_1.1.0 zoo_1.8-8 haven_2.3.1
## [76] ggrepel_0.8.2 cluster_2.1.0 fs_1.5.0
## [79] magrittr_2.0.1 data.table_1.13.2 DO.db_2.9
## [82] lmtest_0.9-38 reprex_0.3.0 RANN_2.6.1
## [85] fitdistrplus_1.1-1 hms_0.5.3 mime_0.9
## [88] evaluate_0.14 xtable_1.8-4 readxl_1.3.1
## [91] gridExtra_2.3 compiler_4.0.3 KernSmooth_2.23-18
## [94] crayon_1.3.4 shadowtext_0.0.7 htmltools_0.5.2
## [97] mgcv_1.8-33 later_1.1.0.1 speedglm_0.3-2
## [100] lubridate_1.7.9.2 DBI_1.1.0 tweenr_1.0.1
## [103] formatR_1.7 dbplyr_2.0.0 MASS_7.3-53
## [106] Matrix_1.2-18 cli_2.2.0 igraph_1.2.6
## [109] pkgconfig_2.0.3 rvcheck_0.1.8 plotly_4.9.2.1
## [112] xml2_1.3.2 bslib_0.3.1 XVector_0.30.0
## [115] rvest_0.3.6 digest_0.6.27 sctransform_0.3.1
## [118] RcppAnnoy_0.0.17 spatstat.data_1.5-2 rmarkdown_2.14
## [121] cellranger_1.1.0 leiden_0.3.5 fastmatch_1.1-0
## [124] uwot_0.1.9 shiny_1.5.0 lifecycle_0.2.0
## [127] nlme_3.1-150 jsonlite_1.7.1 fansi_0.4.1
## [130] pillar_1.4.7 lattice_0.20-41 fastmap_1.1.0
## [133] httr_1.4.2 survival_3.2-7 GO.db_3.12.1
## [136] glue_1.4.2 spatstat_1.64-1 png_0.1-7
## [139] bit_4.0.4 ggforce_0.3.2 stringi_1.5.3
## [142] sass_0.4.1 blob_1.2.1 memoise_1.1.0
## [145] irlba_2.3.3 future.apply_1.6.0